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Abstract 

By means of the Direct Simulation Monte Carlo method, the Boltz- 
mann equation is numerically solved for a gas of hard spheres enclosed 
between two parallel plates kept at different temperatures and subject to 
the action of a gravity field normal to the plates. The profiles of pressure, 
density, temperature and heat flux are seen to be quite sensitive to the 
value of the gravity acceleration g. If the gravity field and the heat fiux 
are parallel > 0), the magnitudes of both the temperature gradient and 
the heat flux are smaller than in the opposite case {g <G). When consid- 
ering the actual heat fiux relative to the value predicted by the Fourier 
law, it is seen that, if p > 0, the ratio increases as the reduced local field 
strength increases, while the opposite happens if gr < 0. The simulation 
results are compared with theoretical predictions for Maxwell molecules. 

1 Introduction 

The steady planar Fourier flow is one of the basic nonequilibrium states. It cor- 
responds to a macroscopic system enclosed between two parallel infinite plates 
located at z = and z — L and kept at different temperatures T_ and T+. 
After a certain transient period, the system reaches a steady state characterized 
by a temperature gradient dT/dz along the direction z normal to the plates 
and a constant heat flux q^. In the case of a fluid described by the Navier- 
Stokes equations, the heat flux is just proportional to the temperature gradient 
(Fourier law): 

riT 

= -<nz)) — , (1) 
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where k{T) is the thermal conductivity coefficient, which, in general, depends 
on the local temperature T. In principle, the validity of the linear relation 
is restricted to small gradients, i.e. to £ ^ X, where £ is the characteristic 
hydrodynamic length defined as ^ = T\dT/dz\~^ and A is a microscopic charac- 
teristic length (such as the mean free path in the case of a dilute gas). However, 
computer simulations of both dense fluids and dilute gases |2], ^ , as well as 
kinetic theory treatments ^, ||, H, 0, § show that equation (|l|) is an excellent ap- 
proximation even ii £ ^ X. In the special case of dilute gases, Asmolov et al. |Q| 
found an exact solution of the Boltzmann equation for Maxwell molecules in 
which equation (|l|) is verified for arbitrary values of the thermal gradient. The 
same result is obtained from an exact solution Q of the Bhatnagar-Gross-Krook 
(BGK) kinetic model for any interaction potential. 

The effect of gravity on the heat conduction of dilute gases is usually ne- 
glected. This is because the characteristic distance associated with gravity (i.e. 
the scale height h = v^/g, where vq is the thermal velocity and g is the gravity 
acceleration) is much larger than £ and A under usual laboratory conditions. 
Thus, the Fourier law (|^) still holds if /i ^ £ ^ A. On the other hand, an 
interesting question from a physical point of view is whether or not the heat 
conduction is influenced by a gravity field g along the z-direction when the 
conditions of rarefaction and/or the strength of the field are such that the ratio 
X/h is not negligibly small. In the case of Earth's atmosphere, for instance, h 
varies only within the range of 5-10 km up to an altitude of 100 km 
On the other hand, A increases from lO"'^ cm at the surface to tens of km at 
an altitude of 500 km. Consequently, X/h ^ 10^^^ at the surface but rapidly 
increases with the altitude, being A//i '--^ 1 at the base of the exosphere [ pT| . 

In the study of the gravity effect on the heat conduction the main quantity 
of interest is the reduced thermal conductivity coefficient K*(e, 5*), where 



A 

1 



ainT 



dz 



(3) 



a* = -= (41 
^ ~ h 2kBT/m' ^ ' 

The quantity e is a dimensionless measure of the thermal gradient over the scale 
of the mean free path, while g* is a measure of the strength of the gravity field 
over the same scale. In equation (^) we have identified the thermal velocity as 
vq = {iksT /m)^/"^ ^ where ks is the Boltzmann constant and m is the mass of 
a particle. As a convenient definition of the mean free path we take 
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where p = nksT is the hydrostatic pressure, n being the number density. The 
above definition of A is based on the result for the thermal conductivity coeffi- 
cient K in the BGK model It must be pointed out that all the quantities 
in (|2|-||) are local, i.e. they depend on z. For instance, k* is the ratio between 
the actual value of the heat flux (which is in fact uniform in the steady state) 
and the local value obtained from (||) with the actual local values of the temper- 
ature and its gradient. If the system is large enough (L — oo) and we restrict 
ourselves to the hulk domain (i.e. far from the boundaries), it is expected that 
all the space dependence of k* occurs through e and g* only. 

The problem of elucidating the effect oi g* ^ Q on the effective thermal 
conductivity k* has been recently addressed by us from a theoretical point of 
view |l^ |l^ . In Rcf. the Boltzmann equation for Maxwell molecules 

was solved by a perturbation expansion through order g^ . The result is 

«*(£, 5*) = '^*(e, 0) -f ^eg*-K (^^-f 503.762^ g*VO(<?*'), (6) 

where K*(e,0) = 1. Henceforth we take the convention that g > i) when the 
field g is parallel to the heat flux vector q and g < Q when g is antiparallel to 
q. In the first case, according to equation (^, the heat transport is enhanced 
with respect to the Navier-Stokes prediction (k* > 1), while it is inhibited if 
g < Q. Of course, in the limit of a negligible field [g* — > 0), the validity of the 
Fourier law for arbitrary e Q is recovered (k* ^ 1). A similar analysis has been 
carried out in the context of the BGK model of the Boltzmann equation, also 
for Maxwell molecules [|l^. This allowed us to perform the expansion through 
sixth order in the field, the behaviour of the numerical coefficients indicating 
that the series expansion is only asymptotic. The result to second order is 

«*(e,5*) = «^*(e,0) -t- ^-^eg* + + 9^ + 0{g*^). (7) 




The asymptotic analysis of Ref . |16f| agrees well with a finite-difference numerical 
solution of the BGK equation |17| . Comparison between equations (^) and (Q) 
shows that the BGK model tends to overestimate the influence of the gravity 
field. A much more complex gravity effect appears in the planar Couette flow 
where normal and shear stresses are present in addition to heat transport. 
Thus far, all the previous studies about the influence of gravity on the heat 
conduction have been restricted to Maxwell molecules. In addition, most of them 
have been based on theoretical asymptotic analyses of the Boltzmann or BGK 
equations, the investigation of Ref. JT^ being the only exception. The main 
merit of the Maxwell interaction is that it usually makes the analytic treatment 
of the Boltzmann equation more manageable, but otherwise it is a rather un- 
realistic potential. Structural, thermodynamic and transport properties of real 
fluids are much better captured by the hard-sphere model ||l9| . The aim of this 
paper is to investigate the gravity effect on the planar Fourier flow in the case of 
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a dilute gas of hard spheres by solving numerically the Boltzmann equation by 
means of the DSMC method |2^. As will be seen, the simulation results agree 
qualitatively with the theoretical analyses. The physical problem is stated in 
section ^ where a special attention is paid to the choice of the boundary condi- 
tions. The simulation method is described in section ||. Section ^ presents the 
results. Finally, the conclusions are summarized in section 0. 



2 Planar Fourier flow in the presence of a grav- 
ity fleld 

Let us consider a dilute gas of hard spheres enclosed between two parallel plates 
located at z = and z = L. Both plates are maintained at temperatures T_ 
and T+, respectively. Without loss of generality, we will assume that T+ > T_. 
In addition, a constant gravity field g = ~gz is applied. Figure Q presents 
a schematic illustration of the system geometry. Under these conditions, the 
Boltzmann equation reads |l^ 

d d d '\ f f 

X [/(z,v',t)/(z,v'i,t) 

-/(z,v,i)/(z,vi,0]. (8) 

Here, f{z,v,t) is the one-particle velocity distribution function, a is the diam- 
eter of a sphere, 6 is the Heaviside step function, g = v — vi is the relative 
velocity, ct is a unit vector in the direction of the line joining the centers of 
the two colliding particles, and v' = v — (g • a)d and v'^ = vi -|- (g • d)d are 
post-collisional velocities. In equation we have already assumed that, due 
to the geometry of the problem, the only relevant space coordinate is z. The 
hydrodynamic fields and the associated fluxes can be expressed as moments of 
the distribution function: 



n{z,t)^ J dvf{z,v,t), (9) 
P{z,t) j dwwwf{z,w,t), p = nfci3T=iTrP, (10) 

qiz^t)^"^ I d^v'vf{z,^,t). (11) 

In the expressions of the pressure tensor, equation ([l0|), and of the heat flux, 
equation (^l]), we have assumed the absence of a flow velocity field. Multiply- 
ing equation (||) by (vz,v'^) and integrating over velocity space, one gets the 
following steady state balance equations 

^P,, + rung = 0, (12) 



4 



^..=0. (13) 

In order to solve the Boltzmann equation (||) one needs to complement it 
with initial and boundary conditions. Since we will focus on the steady state, the 
particular choice of initial condition is not relevant here. As for the boundary 
conditions, they can be expressed in terms of the kernels K±{v,v') defined as 
follows. When a particle with velocity v' hits the wall aX z = L, the probability 
of being reemitted with a velocity v within the range dv is K+{\r,v')dv, the 
kernel K^(v,v') represents the same but at z = 0. The boundary conditions 
are then 

Qi±v,)\v,\f{z^{0,L},^,t) = Q{±v,) j d^r'\v',\K^{^y) 

xe(T«;)/(^ = {0,L},v',t). (14) 

In this paper we consider boundary conditions of complete accommodation with 
the walls, so that K±{'v, v') = K±{y) does not depend on the incoming velocity 
v' and can be written as 

X±(v) = A±le(T^'.)|«.|</'±(v), A± = j dve(Tt'.)|«.|0±(v), (15) 

where (pz^ (v) represents the probability distribution of a fictitious gas in contact 
with the system at z = {0, L}. Equation ( [l5| ) can then be interpreted as meaning 
that when a particle hits a wall, it is absorbed and then replaced by a particle 
leaving the fictitious bath. Of course, any choice of (j>±{v) must be consistent 
with the imposed wall temperatures, i.e. 

kBT± = ^m j d^fv'^4>±{w). (16) 

The simplest and most common choice is that of a Maxwell-Boltzmann distri- 
bution: 

3/2 



These boundary conditions have been frequently used in molecular dynamics 
simulations jl], |^ as well as in kinetic theory analyses ||, |lj, ^ . Under these 
conditions, the system is understood to be enclosed between two independent 
baths at equilibrium at temperatures and , respectively. While the condi- 
tions (|lj) are adequate for analysing boundary effects ||lj, ^ , they are not very 
efficient when the interest lies on the transport properties in the bulk. In order 
to inhibit the influence of boundary effects, it is more convenient to imagine that 
the two fictitious baths are in nonequilibrium states resembling the state of the 
actual gas near the walls. More specifically, we can assume that the fictitious 
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gases are described by the BGK equation, whose exact solution for the steady 
planar Fourier flow (in the absence of gravity) is known M. In that case, 



(v) = 



3/2 



exp 




2TrkBT±J '\ 2kBT± 

/ rfTe((l - r)sgni;;,)r-3/2 
Jo 



{2kBT±/m 

<^±\Vz\ 



,1/2 



X exp 



-i2kBT±/m) 



1/2 



1 



e±Vz 2kBT±T 



(18) 



where we have additionally assumed statistical independence among the three 
velocity components. In equation ( p^ ) Cip is the local reduced thermal gradient 
at z = {0,L}. If we formally take the limit e± 0, equation ( p^ ) reduces to 
equation ( p7| ) On the other hand, if T_ ^ r+, then e± ^ 0. The exact 
solution of the BGK model ||, ^ has the properties dT^^^/dz — const (for hard 
spheres) and p = const; from them, it is easy to obtain 



15 T 



1/2 



T 



1/2 



-1/2 



where 



dz n{z) 



(19) 



(20) 



is the average density. This second class of boundary conditions were first 
proposed in Ref. , where they proved to be much more efficient than the con- 
ditions (17) for the heat conduction problem. Following the same terminology 
as in Ref. [ || , we will refer to the "equilibrium" conditions ( p^ ) as conditions of 
Type I and to the "nonequilibrium" conditions (nsl) as conditions of Type II. 



3 Simulation method 

In order to solve numerically the Boltzmann equation (|^) with both types of 
boundary conditions, we have used the so-called Direct Simulation Monte Carlo 
(DSMC) method (20[ |2^. Comparison with known exact solutions of the Boltz- 
mann equation under strong nonequilibrium conditions [p^ proves the reliability 
and efficiency of the DSMC method to solve the Boltzmann equation. In this 
method, the velocity distribution function is represented by the velocities {v^} 
and positions {zi} of a sufficiently large number of particles N. Given the ge- 
ometry of our problem, the physical system is split into layers of width Az, 
sufficiently smaller than the mean free path. The velocities and coordinates are 
updated from time t to time t + At, where the timestep At is much smaller than 
the mean free time, by applying a convection step followed by a collision step. 
In the convection step, the particles are moved ballistically, i.e. Viz Viz — gAt 
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and Zi ^ Zi + VizAt — ^g{At)^. In addition, those particles crossing the bound- 
aries are reentered with velocities sampled from the corresponding probability 
distribution K±{v). The collision step proceeds as follows [^0|, For each 
layer a, a pair of potential collision partners, i and j, and a unit vector aij 
are chosen at random with equiprobability. The collision between particles i 
and j is then accepted with a probability equal to Q{gij ■ aij)ujij /uini&K, where 
gij = Vi — Vj is the relative velocity, ujij = (g^ • (Ty)47rcr^rtQ, is the collision rate 
of the pair {i,j), Ua being the number density in layer a, and o^max is an upper 
bound estimate of the collision rate in the layer. If the collision is accepted, 
post-collisional velocities v- ,^ = v^j =F {gij ■ ai)ai are assigned to both parti- 
cles. After the collision is processed or if the pair is rejected, the routine moves 
again to the choice of a new pair until the required number of candidate pairs 
iiVcWiiiaxAi in the layer, where Na is the total number of particles in layer a, 
has been processed. 

In the course of the simulations, the following "coarse-grained" local quan- 
tities are computed. The number density in layer a is 

- ^e.(.4 (21) 



{N/L)Az NAz ^ 

where Qa{z) is the characteristic function of layer a, i.e. Qa{z) — 1 if z belongs 
to layer a and is zero otherwise. Similarly, the temperature, the pressure tensor 
and the heat flux of layer a are 

N 

fci3T„ = ^ = -^5:e„(z,).f, (22) 

ria 3iVa ^ 

2—1 



L ^ 

-^^m^eo(20v4V„ (23) 

4=1 

T ^ 

'"^e„(zO^;?v4. (24) 



NAz 2 

4=1 



According to equation ([l^), q is a constant in the steady state. Thus, wc have 
also evaluated the average heat flux as 



A 1 ^ 

Az \ ^ 1 TO 



L N 2 

a 4=1 



The standard definition of mean free path in the case of hard spheres is [ p^ 

\/2mTa'^ ^ ^ 
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Quantity 


Unit 


Reference value 


Temperature (T) 
Mass (m) 


r+ 
m 


10-' K 

3 X 10-26 kg 


Length (z, L) 


x' 


625 m 


Speed [v) 


(2fcBr+/m)i/2 
A'(2fcBT+/m)-i/2 

A'"'(2fcsT+/m) 


959 ms-^ 


Time (t) 


0.652 s 


Acceleration (g) 


1.47 X 10^ ms-2 


Number density (n) 


n 


4 X 10^^ m-3 


Pressure {p, Pzz) 


n{2kBT+) 


1.10 X 10-'' Nm-2 


Heat flux (g^) 




0.106 Jm-2s-i 



Table 1: Units used in this paper for the relevant quantities. The third column 
gives the "real" values of those units taking (t = 3A, m = 3x 10-26 
T+ = 10^ K and n = 4 x 10^^ m""^ as a reference example. 



This quantity is not exactly the same as the mean free path defined by equation 
which is based on the BGK model. The thermal conductivity of a dilute 
system of hard spheres is (in the first approximation) p3| 

_ 75(7rmfcBr)i/2fcg 

so that A — (15-\/7r/16)A'. In the following, we take units such that A = 
(•\/2n7rcr2)-^ = 1 (length unit), to = 1 (mass unit), (2/cbT+/to)^/2 — i (speed 
unit), T_|_ = 1 (temperature unit) and n — 1 (density level). The units of 
these and other related quantities are given in Table ^ The table also gives 
the values of those units taking as a reference example a gas with cr = 3 A, 
TO = 3 x 10-26 kg, T+ = 10^ K and n = 4 X 10^^ m'^, which are typical of the 
atmospheric conditions at an altitude of 220 km p^ . 

Each physical situation is then characterized by the values of L, T_ and 
g only. In terms of quantities expressed in the above units, the reduced local 
parameters e and g* , defined by equations (||) and (^, are given as 



16pTi/2 



dz 



(28) 



9* = (29) 
^ 32 p ^ ' 



According to the Fourier law, equation (hi), the heat flux is 



'^^ " 64 5z ' ^"^"^ 

where we have taken into account that k oc T^/2 fQj- ^ hard-sphere gas and have 
used the identity T^/^dT = ^dT^^^. Equation (|^ imphes that in the steady 
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state described by the Fourier law one would have dT^^'^/dz — const. While 
this is not in general true, the slope of T^/^ is expected to change more smoothly 
than that of T. This is why in equation ( p8| ) e is written in terms of the former 
rather than in terms of the latter. 

In the simulations we have taken N 5000 particles, a layer width Az = 0.1 
and a timestep At = 0.004. We have started from initial conditions of the form 



/(z,v,0) =n(z,0) 



27rfcBr(z,0) 



3/2 



exp 



2fcBT(z,0) 



(31) 



where n(z,0) cx 1/T(z,0) andr(z,0) = T_(l+cz/L)2/3, with c = {T+/T^f/^~ 
1. After a time period t = 200 the system has already relaxed to the steady 
state We follow the evolution of the system until t — 2000 and average 
the relevant quantities over 5000 snapshots equally spaced between t = 200 and 
t = 2000. 



4 Results 

By using the method outlined in the previous section, we have analysed 51 
different states. In all of them, the separation between the plates has been taken 
as L — 10. Two different temperature ratios have been considered {T^ = 0.01 
and T_ = 0.05) and boundary conditions of Types I and II have been applied. 
For each one of these four combinations, 12 or 13 different values of g have been 
taken, typically in the range —0.024 < g < 0.014. We illustrate in figures ||-|| the 
profiles found in the simulations by choosing the case T_ = 0.05 with boundary 
conditions of Type II as a reference example. The corresponding temperature 
profiles are shown in figure | for 5 = -0.016, -0.008, 0, 0.008 and 0.012. We 
observe that the temperature gradient is larger for g < (i.e. when the gravity 
field is antiparallel to the heat flux) than for g > (gravity field parallel to 
the heat flux). This implies that the magnitude of the heat flux is expected to 
be larger in the former case than in the latter. This is confirmed by figure ^, 
where the profile of is plotted for the same situations as in figure |^. Figure 
^ also shows that, except for statistical fiuctuations, the results are consistent 
with Qz = const. This is a consistency test that a steady state has been reached 
in the simulations [cf. equation (p^j. The fact that the heat fluxes are constant 
and yet the profiles of T^^-^ are nonlinear can be traced back to local deviations 
from the Fourier law. 

The profiles of p and Pzz are shown in figure ^. At g = the simulation 
results are consistent with a constant Pzz- On the other hand, in agreement 
with equation (|l2|) , Pzz is an increasing (decreasing) function of z when g < 
{g > 0). The hydrostatic pressure p is slightly larger than Pzz for z/L smaller 
than about 0.3-0.4, while it is slightly smaller than Pzz for larger distances 
from the cold wall. It is worthwhile noting that the local density n = 2p/T is 
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much smaller in the hotter layers than in the colder ones, as seen in figure 
Nevertheless, this disparity in the population of particles is widely influenced 
by the sign of g. Thus, for g = —0.016 the densities near the cold and hot walls 
are n ~ 4 and n ~ 0.6, respectively, while those values are n ~ 8 and n ~ 0.3 for 
g = 0.012. The large densities near the cold wall are responsible for the abrupt 
change of the pressure in that region, in agreement with the balance equation 

In the cases T_ = 0.05 (with boundary conditions of Type I) and T_ = 0.01 
(with boundary conditions of Types I and II) we have obtained results similar 
to those displayed in figures The most interesting difference is that, as 

expected, the boundary effects are much less important when applying boundary 
conditions of Type II than those of Type I. This is illustrated in figure ^, where 
it can be seen that, given a value of L, T_ and g, the jump temperature at 
the walls is much smaller in the case of boundary conditions of Type II. This 
extends to g 7^ the observation made in Ref. Q for 5 = 0. 

As stated in section 0, we are mainly interested in investigating the influence 
of the gravity strength on the heat flux relative to the value predicted by the 
Navier-Stokes approximation with the actual thermal gradient. This ratio is the 
reduced thermal conductivity k*, equation (|^). In the limit in which boundary 
effects are negligible {L 00), it is expected that k* depends on position 
only through a functional dependence on the local parameters e and g*, defined 
by equations ^ and (^) or, equivalently, by equations (|2^) and ( p9| ) in our 
units. In order to minimize as much as possible the unavoidable boundary 
effects associated with a finite L, we will focus on the region around a point zq 
sufficiently far from both boundaries. By measuring the pressure, temperature 
and thermal gradient at z — zq (the heat flux is measured as the average value 
g^), we can compute the associated values of k*, g* and e = ep. The question 
arises as to how to choose the value of zq. Here we have applied the following 
criterion. First, the value of zo for g = is such that the expected number of 
collisions a particle experiences when going from z = to z = zq or from z ~ L 
to z = Zq is 5, i.e. Af{zQ) — 5, where 

r dz' , , 

Notice that, since A' ~ 1/n in our units, M(L) = L = 10. In this way we 
have obtained eo = 0.33 (T_ = 0.05, Type I), cq = 0.36 (r_ = 0.05, Type 
II), eo = 0.44 (r_ = 0.01, Type I) and eo = 0.48 (T_ = 0.05, Type II), aU of 
them with g — If we followed the same method to choose zo when g 0, 
then we would obtain a different value of eo each time and that would hinder 
the analysis about the direct influence of the gravity field on the coefficient k* . 
Therefore, as the second part of the criterion, we fix the above values of eo 
for each one of the four combinations of r_ and type of boundary conditions 
and then determine zq to accommodate the corresponding eo. The values of zo, 
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Q 


zqI L 


7V(zo) 






K* 


Ak* 


-0.024 


2.55 


3.3 


-0.095 


0, 


.903 


-0.040 


-0.020 


2.55 


3.5 


-0.079 


0, 


.905 


-0.038 


-0.016 


2.60 


3.8 


-0.063 


0, 


.906 


-0.037 


-0.012 


2.60 


4.1 


-0.048 


0, 


.917 


-0.026 


-0.010 


2.60 


4.3 


-0.040 


0, 


.921 


-0.022 


-0.008 


2.55 


4.3 


-0.032 


0, 


.927 


-0.016 


-0.006 


2.60 


4.5 


-0.025 


0. 


.934 


-0.009 


-0.004 


2.50 


4.5 


-0.017 


0, 


.939 


-0.004 


0.000 


2.70 


5.0 


0.000 


0, 


.943 


0.000 


0.003 


2.80 


5.3 


0.013 


0, 


.951 


0.008 


0.006 


2.85 


5.6 


0.028 


0, 


.969 


0.026 


0.008 


2.90 


5.8 


0.039 


0, 


.978 


0.035 


0.010 


3.20 


6.3 


0.052 


0, 


.986 


0.043 



Table 2: Values of zq, M(z{^\ g* , k* and An* for different values of g in the 
case r_ = 0.05 with boundary conditions of Type I. The values of zq are such 
that the reduced thermal gradient is eg = 0.33. 

Af{zo), g* and k* obtained in this way are given in Tables |-| As we can see, 
the point zq is always closer to the cold wall than to the hot wall. However, 
when the separation is measured in terms of the number of collisions rather than 
as an actual distance, it turns out that the point zq is "closer" to the hot wall, 
i.e. Af{zo) > 5, if g > 0, while the opposite happens if g < 0. 

In the absence of gravity {g = 0), we observe that k* is smaller than 1. 
This is in part due to a residual influence of boundary effects , as indicated 
by the fact that k* is closer to 1 with boundary conditions of Type II than 
with those of Type I. From the exact solution of the Boltzmann equation for an 
unbounded system of Maxwell molecules [Q, it follows that k* — 1 for g* = 
and arbitrary e. However, it is possible that K*(e,0) is slightly smaller than 
1 in the case of hard spheres, even if any boundary effect is removed. Thus, 
in order to characterize the pure gravity-dependence of the effective thermal 
conductivity k* at a given value of e, we define 

A^*(6,g*)^ 5*) -«*(6,0). (33) 

The values of Ak* for the corresponding fixed values of cq are also included in 
Tables ^-|^ and are plotted in figure 0. It is observed in figure that, although 
somewhat scattered, the points tend to lie on smooth curves. The behaviour 
of Ak* is in qualitative agreement with the theoretical predictions for Maxwell 
molecules, equations (^) and (|^). More specifically, the departure from the 
Fourier law, as measured by Ak*, is positive when the gravity field is parallel 
to the heat flux [g* > 0), while it is negative when both vectors are mutually 
antiparallel (g* < 0). In addition, the magnitude of the deviation is larger 
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.9 


zo 


,/L 


AA(zo) 


.9* 






Ak* 


-0.020 


2, 


.45 


3.9 


-0.068 





.908 


-0.055 


-0.016 


2, 


.45 


4.1 


-0.055 





.921 


-0.042 


-0.012 


2, 


.45 


4.3 


-0.042 





.934 


-0.029 


-0.010 


2, 


.40 


4.4 


-0.035 





.940 


-0.023 


-0.006 


2, 


.50 


4.7 


-0.022 





.949 


-0.014 


-0.003 


2, 


.50 


4.9 


-0.011 





.955 


-0.008 


0.000 


2, 


.50 


5.0 


0.000 





.963 


0.000 


0.004 


2, 


.60 


5.5 


0.016 





.983 


0.020 


0.008 


2, 


.60 


5.5 


0.032 





.988 


0.025 


0.010 


2, 


.90 


6.2 


0.046 


1 


.016 


0.053 


0.012 


3, 


.00 


6.4 


0.058 


1 


.036 


0.073 


0.014 


3, 


.55 


7.0 


0.077 


1 


.064 


0.101 



Table 3: Values of zqj A/'(zo), 5*, n* and Ak* for different values of g in tlie 
case r_ = 0.05 with boundary conditions of Type 11. The values of zq are such 
that the reduced thermal gradient is eg = 0.36. 
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zo/L 


N{ 


z,) 


9* 




K* 


Ak* 


-0.024 


2.00 


3, 


,1 


-0.108 


0, 


.833 


-0.066 


-0.020 


2.00 


3, 


.4 


-0.091 


0. 


.838 


-0.061 


-0.016 


2.00 


3, 


.6 


-0.073 


0, 


.838 


-0.061 


-0.014 


2.00 


3, 


.8 


-0.065 


0, 


.849 


-0.050 


-0.012 


1.95 


3. 


.9 


-0.056 


0. 


.861 


-0.038 


-0.008 


2.00 


4, 


.3 


-0.038 


0, 


.871 


-0.028 


-0.004 


2.00 


4, 


.6 


-0.020 


0, 


.880 


-0.019 


0.000 


2.05 


5. 


.0 


0.000 


0. 


.899 


0.000 


0.004 


2.15 


5. 


.5 


0.022 


0, 


.917 


0.018 


0.006 


2.40 


6, 


.0 


0.036 


0, 


.921 


0.022 


0.008 


2.70 


6. 


.5 


0.052 


0. 


.930 


0.031 


0.010 


3.00 


7. 


.0 


0.074 


0. 


.964 


0.065 


0.011 


3.50 


7. 


.5 


0.092 


0. 


.994 


0.095 



Table 4: Values of zq, Af{zo), g*, k* and Ak* for different values of g in the 
case r_ = 0.01 with boundary conditions of Type I. The values of zq are such 
that the reduced thermal gradient is eo = 0.44. 
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Q 


zqI L 


7V(zo) 




K* 


Ak* 


-0.024 


1.85 


3.4 


-0.088 


0.837 


-0.075 


-0.020 


1.80 


3.6 


-0.075 


0.851 


-0.061 


-0.016 


1.80 


3.8 


-0.061 


0.864 


-0.048 


-0.012 


1.80 


4.1 


-0.047 


0.879 


-0.033 


-0.008 


1.80 


4.4 


-0.032 


0.888 


-0.024 


-0.004 


1.85 


4.7 


-0.016 


0.890 


-0.022 


0.000 


1.90 


5.0 


0.000 


0.912 


0.000 


0.004 


1.95 


5.5 


0.019 


0.939 


0.027 


0.006 


2.00 


5.7 


0.030 


0.943 


0.031 


0.008 


2.10 


6.0 


0.043 


0.974 


0.062 


0.010 


2.20 


6.3 


0.058 


0.992 


0.080 


0.011 


2.30 


6.5 


0.068 


1.010 


0.098 


0.012 


2.70 


7.1 


0.084 


1.047 


0.135 



Table 5: Values of zq, A/'(2:o), 5*, '«* and Ak* for different values of g in the 
case r_ = 0.01 with boundary conditions of Type II. The values of zq are such 
that the reduced thermal gradient is eg = 0.48. 



in the former case than in the latter, i.e. AK*(e,(7*) > — AK*(e, — g*). From 
figure we can also conclude that the gravity influence is less dramatic when 
the boundary effects are more important, since |Ak*| tends to be smaller in 
the case of boundary conditions of Type I. From the data corresponding to the 
boundary conditions of Type II we can estimate that AK*(e, y*) ~ Beg* + ■ • •, 
where the coefficient B has a value between 2 and 2.5. This is about 4 times 
smaller than the exact coefficient B = 46/5 obtained in the case of Maxwell 
molecules. We are not able at present to elucidate which part of this discrepancy 
in the numerical value of B is attributable to boundary effects still present in our 
simulation results and which part is due to the role played by the interaction. 
Besides, when comparing results obtained with different interactions, one must 
have in mind that the choice of appropriate dimcnsionless parameters is not 
unique. In our case, we have defined e, equation (|^), and g*, equation by 
using the mean free path given by equation (|^), which is based on the BGK 
model. If, on the other hand, we had used the standard mean free path of hard 
spheres, equation (26), then we would have AK*{e,g*) ~ B'eg* + ■ • ■, where 
B' = {X/X'fB ~ 2. 765. 



5 Conclusions 

In this paper we have numerically solved the Boltzmann equation (by means of 
the DSMC method) for a steady heat conduction problem of hard spheres in 
the presence of a gravity field. The gas is enclosed between two parallel plates 
separated a distance equal to L = 10 (average) mean free paths. Two different 
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temperature ratios {T-/T+ = 0.01 and T^/T+ = 0.05) and two alternative types 
of boundary conditions (boundary effects being more important in the case of 
Type I than in the case of Type II) have been considered. For each one of these 
four possibihties we have apphed 12 or 13 different values of a constant gravity 
field g = —gz normal to the plates. The sign criterion is such that .g > means 
a field antiparallel to the thermal gradient (and hence parallel to the heat flux 
vector), while g < means the opposite. 

The first conclusion we draw from the results is that the hydrodynamic 
profiles are rather sensitive to the value of g. While the pressure p is roughly 
uniform if 17 = 0, it decreases (increases) with z ii g > {g < 0), this effect 
being more important as the magnitude of g grows. This is not surprising since 
it is an extension to nonequilibrium states of the equilibrium barometric law 
p{z) (X Q-'^gz/kBT ^ ^ jggg obvious effect appears in the case of the temperature 
profile. If there were no temperature jumps at the walls, the temperature of 
the gas would change from T_ at 2 = to T_|_ at z = L, irrespective of the 
value of g. However, due to unavoidable boundary effects, r(0) > T„ and 
T{L) < r+. Our simulation results show that the temperature jump at the cold 
(hot) wall decreases as the value of 5 increases (decreases). For instance, in the 
case r_ = 0.05 with boundary conditions of Type 11, r(0) - T_ ~ 0.04 and 
T+ - T{L) ~ 0.04 for g ^ -0.016, while T(0) - T_ ~ 0.01 and T+ - T{L) ~ 0.20 
for g = 0.012. In other words, the temperature jump at a wall decreases as the 
relative density near that wall increases. The more positive (negative) g is, the 
larger the density is near the cold (hot) wall and the smaller the temperature 
jump is. Since the temperature jump is more important near the hot wall, a 
side effect of the above discussion is that the (average) thermal gradient across 
the system increases when g decreases, so that it is larger for g > than for 
<; < 0. A larger thermal gradient implies a larger magnitude of the heat flux 
and this expectation is confirmed by our simulation results, which show that 
\qz\ clearly increases as g decreases. 

Our interest has not focused, however, on the absolute change of the heat 
flux due to the gravity field, but on its change relative to the Navier-Stokes 
prediction when the actual temperature gradient is considered. To that end we 
have introduced the ratio k* defined by equation (^), which is a local quantity 
that in the bulk region is expected to depend on space only through a func- 
tional dependence on the reduced thermal gradient e, equation (^, and gravity 
strength g* , equation (^). For each one of the four different combinations of 
T^/T-^- and boundary conditions we have fixed a value e = eo (eg ~ 0.3-0.5) 
and have analysed the (7*-dependence of AK*{eo,g*) = K*{eo,g*) — K*(eo,0). 
The simulation results for hard spheres presented in this paper are in qualita- 
tive agreement with those obtained for Maxwell molecules by a perturbation 
analysis of the Boltzmann |l^ and the BGK ||l^ equations. More specifically, 
when the field and the heat flux are parallel {g* > 0) the gravity induces an 
enhancement of the (relative) heat conduction {An* > 0), while the opposite 
happens when both vectors are mutually antiparallel. In addition, the influence 
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of gravity is more pronounced in the former case than in the latter. At a quan- 
titative level, on the other hand, the values of |A«;*| reported here are typically 
smaller than those theoretically estimated for comparable values of e and g*. A 
certain part of this difference is possibly due to boundary effects still present 
in our simulations and absent in the theoretical analyses of Refs. and p6{ . 
This expectation is supported by the fact that |Ak*| is generally smaller in the 
case of conditions of Type I than in that of Type II, thus indicating that bound- 
ary effects tend to mitigate the influence of gravity. Notwithstanding this, the 
remaining difference leads us to conclude that the response of the system to the 
application of the field, as measured by Ak*, is less important in the case of 
hard spheres than in the case of Maxwell molecules. 
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Figure 1: Schematic illustration of the system geometry. 



Figure 2: Profile of T^/^ in the case T_ = 0.05 with boundary conditions of 
Type II. The values of g are, from top to bottom, g = -0.016, -0.008, 0, 0.008 
and 0.012. The labels C and H denote the locations of the cold and hot walls, 
respectively. 



Figure 3: Profile of in the case T_ = 0.05 with bomidary conditions of Type 
II. The values of g are, from bottom to top, g = -0.016, -0.008, 0, 0.008 
and 0.012. The labels C and H denote the locations of the cold and hot walls, 
respectively. 



Figure 4: Profiles of p (solid lines) and Pzz (dashed lines) in the case T_ = 0.05 
with boundary conditions of Type II. The values of g arc, from top to bottom 
at the right end, g = -0.016, -0.008, 0, 0.008 and 0.012. The labels C and H 
denote the locations of the cold and hot walls, respectively. 



Figure 5: Profile of n in the case T_ = 0.05 with boundary conditions of Type 
II. The values of g are, from top to bottom at the right end, g = —0.016, —0.008, 
0, 0.008 and 0.012. Note that the vertical axis is in logarithmic scale. The labels 
C and H denote the locations of the cold and hot walls, respectively. 



Figure 6: Profile of T^/^ in the case T_ = 0.01 with boundary conditions of 
Types I (dashed lines) and II (solid lines). The values of g are (a) g = —0.016 
and (b) g = 0.012. The labels C and H denote the locations of the cold and hot 
walls, respectively. 



Figure 7: Plot of Ak* versus the reduced field strength g* for T_ = 0.05 with 
boundary conditions of Type I (open squares, e ~ 0.33), T_ = 0.05 with bound- 
ary conditions of Type II (solid squares, e = 0.36), T_ = 0.01 with boundary 
conditions of Type I (open circles, e = 0.44) and T_ = 0.01 with boundary 
conditions of Type II (solid circles, e = 0.48). The lines are polynomial fits to 
guide the eye. 
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